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This paper presents a numerical analysis of the atmospheric air transonic flow through de Laval nozzles. By na- 


ture, atmospheric air always contains a certain amount of water vapor. The calculations were made using a Laval 


nozzle with a high expansion rate and a convergent-divergent (CD) “half-nozzle”, referred to as a transonic dif- 


fuser, with a much slower expansion rate. The calculations were performed using an in-house CFD code. The 


computational model made it possible to simulate the formation of the liquid phase due to spontaneous condensa- 


tion of water vapor contained in moist air. The transonic flow calculations also take account of the presence of a 


normal shock wave in the nozzle supersonic part to analyze the effect of the liquid phase evaporation. 
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Introduction 


Moist air is a solution (or mixture) of dry air and su- 
perheated water vapor or dry saturated vapor, or dry sa- 
turated water vapor and liquid or ice mist. Depending on 
the combination, three cases can be distinguished. The 
first one is moist non-saturated air, the second — moist 
saturated air and the third — moist oversaturated air. The 
amount of superheated and saturated water vapor in air is 
limited. The vapor contained in air starts to condense if 
the air temperature is lowered to the vapor saturation 
temperature for its component pressure in moist air, 
which is referred to as the dew point. When condensation 
occurs in the moist air transonic flow field, the flow is 
affected by latent heat released in the process. Within this 
study, a condensing flow was produced by moist air ex- 
pansion in nozzles. This was accompanied by formation 
of shock waves in the supersonic parts of the flow field 
due to an increase in outlet pressure. The presented nu- 
merical investigations were carried out to show the effect 
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of on-shock wave evaporation of the liquid phase result- 
ing from water vapor condensation. 

The numerical analyses presented herein were con- 
ducted using convergent-divergent nozzles. A Laval “full 
nozzle” and a “half-nozzle” were used. In the Laval 
“half-nozzle” the upper wall is shaped like a typical con- 
vergent-divergent nozzle, whereas the bottom wall is flat 
along its entire length. The flow through a convergent- 
divergent nozzle is one of the fundamental and most 
common transonic flows. Laval nozzles are widely used 
in technology, e.g. in transport and power engineering, 
mainly to accelerate gas to supersonic speeds. Due to very 
good insulating properties of transonic flows, de Laval 
nozzles are also used in high-voltage circuit breakers. 

Air moisture has an essential impact on many pro- 
cesses taking place in power engineering equipment and 
machinery. If air cannot be dried at the nozzle inlet in 
technological applications using convergent-divergent 
nozzles, it is absolutely necessary to take account of the 
liquid phase formation, which occurs due to condensation 
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of water vapor contained in atmospheric air. During the 
flow through CD nozzles normal shock waves are gener- 
ated due to frequent changes in parameters at the nozzle 
outlet. On these waves, which are mainly normal, the 
flow parameters change dramatically. 

The problem of the moist air condensing flow model- 
ing has been investigated by numerous researchers for 
many years. Both experimental and numerical studies 
have been made in this field. The most known of them 
are the works of Schnerr and his research team [12, 13, 
14]. The works focus both on the analysis of moist air 
internal flows [1, 2] as well as on aerodynamic applica- 
tions such as the flow around an aircraft wing profile [3]. 
Recent works in this field concern external flows in the 
first place and are related to aerodynamic issues [9]. 

The calculations presented herein were made using an 
in-house CFD code, which has been developed and used 
in many engineering applications for over 15 years in the 
moist air and steam condensing flow calculations, to pre- 
dict condensation on the flow fields with shock waves, 
where the Navier-Stokes equations are solved numeri- 
cally using a 3rd-order MUSCL type TVD finite-differ- 
ence scheme with the 2nd-order Runga-Kutta method for 
time integration. The two-equation k-w SST turbulence 
model was selected for the purpose of turbulence modeling. 


Moist Air 


Gases, air included, hardly ever occur in nature in the 
dry state. What is referred to as moist air is a mixture of 
dry air with water vapour or condensed water. The mass 
of moist air is defined as the sum of dry air mass, m,, and 
the total mass of water vapour and water, m, o: 

m=m,+m,)=m, +m, +m, (1) 

The state of moist air (moist gas) is defined by stating 
its moisture (relative or absolute). 

Absolute moisture is defined as the ratio of the water 
vapour mass to the dry air mass: 


; (2) 


a 


x= 
m 
Relative moisture is referred to as the real-to-maxi- 
mum moisture content ratio; it is defined as the ratio of 
the water vapour pressure to the water vapour saturation 
pressure. However, the relative moisture value is often 
given for integral parameters (in the computational prac- 
tice — for stationary parameters at the inlet) and it is ex- 
pressed in percentages. The result is then as follows: 


Ga Py,0 
= 
P VARY (T 0) 
Knowing the relative moisture value and the moist air 


pressure and total temperature, it is possible to find the 
absolute moisture value from the following equation: 


-100% (3) 
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In addition to the moist air relative and absolute 
moisture values, the wetness mass fraction, i.e. the ratio 
of the water mass (condensate mass) to the moist air 
mass, is also defined: 


yoot (5) 
m 
The wetness mass fraction reaches its maximum value 
upon complete condensation of the entire water vapour 
contained in air. 


M, 0 


(6) 


The relationship between absolute moisture and the 
maximum wetness mass fraction is expressed in the fol- 
lowing way: 


Y max = 


mM, o M, o k m Ymax 

Le ____/m___*mx_ (7) 

Ma m— My, o {= mM, o 1- Y max 
m 

In aerodynamic issues air may with a good approxi- 
mation be treated as a perfect gas because under rela- 
tively small pressures and in higher temperatures it be- 
haves like perfect gases. 

In the case of moist air, it may with a close approxi- 
mation be assumed that the temperatures of the mixture, 
the air and the water vapour are equal: 

T=7,=T, (8) 

whereas the pressure of the mixture, according to 

Dalton's law, is equal to the sum of the pressures of dry 
air and water vapour: 


x= 


P= Pat Py (9) 

where, assuming that moist air is described by the 

perfect gas equation of state and the water vapour tem- 

perature is equal to the temperature of dry air, the water 

vapour pressure can be determined from the following 
equation: 


(Ymy) 


Za .(1— Ymax )+ Ymax- y) 


The water vapour and the dry air density may be found 
using the following relations: 
Py = (Ymax -1) p 
Pa =Z (l= ymax) P 
whereas the air and water vapour mixture density p, 
ignoring the liquid phase volume, is calculated from the 
following equation: 


Pat Py (10) 
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The other thermodynamic quantities of moist air are 
found using the following relations: 
h =h, (l- Vmax) + Ay -( max —y)+h, Y 
S= S4 (l- Vmax) + s, -( max —y)+s; Y 
The surface tension for water enveloped in air is de- 
termined experimentally. One of the ways to do so is to 
apply the formula put forward by Pruppacher and Klett 
[10], which for a flat surface of the water-air phase sepa- 
ration may be written as: 


(76.1 + 0.155- (273.15 - T) -107° 
for T > 249.39K 


(11) 


o(T )= 12 
7) (1.1313-3.7091-10°3 -7)-74 10-4 ae 
— 5.6464 
for T < 249.39K 


The moist air dynamic viscosity may be described us- 
ing the Sutherland formula for the perfect gas: 


-17110527341104 Ty" (13) 
in T+1104 | 273 


In the case of moist air for the temperature range of 
200K to 300K, latent heat may be expressed as being 
linearly dependent on temperature [1]: 

L(T) = 3105913.99 —2212.97-T (14) 


Homogeneous Condensation Model 


Upstream the condensation zone the mixture is in a 
non-equilibrium, metastable state. Departures from equi- 
librium are represented in terms of either subcooling AT 
or supersaturation ratio S. Subcooling is defined as AT = 
T,(p,) — T, where T,(p,) is the saturation temperature 
corresponding to the water vapour local partial pressure 
p, and T is the actual moist air temperature (for super- 
heated states AT < 0 and for subcooled states AT >0). 
The supersaturation ratio is defined as S = p,/p,(7) (for 
superheated states S < / and for subcooled states S > /). 
The methods of describing metastable states using either 
AT or S are equivalent and their choice depends on per- 
sonal preferences. Downstream the condensation zone, 
where equilibrium is established, AT = 0 and S =1. 

Due to expansion below the saturation line, the water 
vapour contained in moist air is strongly subcooled (47), 
which leads to the formation of accidental clusters of 
water molecules and, later on as subcooling proceeds, to 
a rise in the number and probability of formation of 
critical droplets. Critical droplets create condensation 
nuclei on which water vapour condenses rapidly. This 
phenomenon is referred to as spontaneous (homogene- 
ous) condensation, which results in primary droplets and 
then subcooled water vapour changes into mist. As the 
expansion rate gets higher, the number and average size 
of primary droplets decrease. The higher the expansion 
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rate, the later spontaneous condensation takes place. Wa- 
ter molecules settle on the droplets being formed and the 
droplets get bigger. In the model of the condensation 
process applied in the calculations presented herein it is 
assumed that secondary effects of condensation, i.e. the 
settlement of droplets on solid surfaces limiting the flow, 
are omitted. 

The liquid phase is represented by a countless number 
of droplets for which rotational motion and internal mass 
circulation are ignored; it is also assumed that there is no 
interaction between the droplets. The liquid phase is in- 
compressible and it is made of spherical droplets. 

Another assumption is that the two-phase flow de- 
scribed herein is a non-slip one. Consequently, the flow 
conservation equations and the equation for the two- 
equation viscous turbulence model are written for a mix- 
ture [4, 15]. Additional conservation equations describing 
the liquid phase formation and growth due to homoge- 
neous condensation, assuming an averaged radius value, 
may be reduced to a system of transport equations for 
two additional conservation variables: moisture p'Yrom 
and the number of droplets in a kilogram of moist air 
P` nhom 

Additional transport equations modelling the liquid 
phase formation and growth due to homogeneous con- 
densation have the following form: 


hom , Enom | Chom , ÊGhom , Shom =0 (15) 
ôt Ox oy Oz 
where Q is the vector of conservation variables, £, F 
and G are flux vectors, and S is the source term. For Car- 
tesian coordinates the column vectors in the equation 
above are defined as follows: 


|p ‘Yh a 
Oh ei = ° ? 
[P ‘Ahom | 
[P Fiom 2 | 0| 
Ehom = ° > E; hom -| > 
LP Phom ‘Y 0 | 


P Nom’ V 0 | 

P’ Yhom’ w| 0 | 
Ghom = ° ? Gs hom = | ? 

P Mom * W 0 | 


4 *3 2 dr 
-Senp-ar(r *Jnom + 3°" Mhom hom’ - 


=P: J hom 
Wetness mass fraction y;,,, and the number of droplets 
per mixture mass unit 7), determined from the equations 
make it possible to find the averaged radius of the drop- 


lets: 
3 1/3 
Thom = [3 | (16) 
An Pi" hom 
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The homogenous condensation model requires a de- 
scription of two mechanisms: nucleation (formation of 
condensation nuclei with critical radius r*) and droplets 
growth. 


Nucleation rate 


The fundamental condition that decides whether the 
condensation nuclei formed in water vapour will stand a 
chance to survive and not evaporate results from a com- 
parison of specific free enthalpy of nuclei formation to 
the enthalpy maximum value. The maximum value of 
specific free enthalpy of nuclei formation is characterized 
by an unsteady condition of equilibrium between the liq- 
uid and the gaseous phase. The condition makes it possi- 
ble to determine the value of the critical radius that the 
arising nuclei should have to survive. 


r° > d8(D.T)y|p cone, = (Ps) 
where g(p,,v,)|, =v,dp. 
Knowing that the pressure on the liquid convex sur- 


face is bigger than on the flat one, according to the fol- 
lowing relation: 


(17) 


| =const 


Pi = Pyt Po =p 


specific free enthalpy for the liquid may be written as: 


2-o(T,) 
g(Pı: Y1) = Yıdpı = vıdp, + va 220) 
Comparing it to free enthalpy for water vapour, the 
following is obtained: 


2-o0(T. 
v,dp,, = vd, + a 2) 
r 


and finally, after transformations: 


v,dp, —vidp, = vd 229) i 
r 

Substituting v, determined from the gas equation of 
state and integrating both sides — the left one from p, to 
p, and the right one from r—»œ tor’, the result is the 
Kelvin equation for the droplet critical radius, i.e. for a 
radius for which the droplet can grow further. Assuming 
the perfect gas model, the equation takes the form 
well-known from literature: 


m 2-0 (18) 


pı: R-T, nf] 
Ps 


Based on the assumed condition of equilibrium be- 
tween the nuclei and the surrounding gaseous phase, the 
condensation nuclei formation rate, i.e. the nucleation 
rate, may be determined in a unit of time and in a unit of 
mass of the mixture J;,.,,. Using the classic theory of nu- 
cleation [5, 6, 8], the nucleation rate may be calculated 
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from the following equation: 
2. 
Jion =e: a m” 
m 
2 | (19) 
. r. . O 


PI 3-k-T, v 

The nucleation rate is very sensitive to changes in 
subcooling (AT) or supersaturation (S). As these quanti- 
ties decide about the critical radius, which in Equation 
(19) is in the second power, even a slight change in sub- 
cooling involves great changes in the rate of critical 
droplets formation. 

In fact there is a difference between the temperature of 
the droplet being formed and the temperature of the 
gaseous phase surrounding it, which has a direct impact 
on the value of the nucleation rate. The effect may be 
taken into account using what is referred to as the non- 
isothermal coefficient put forward by Kantrowitz [7]: 


-1 
ae] oy 
y+1 R-T\R-T, 2 

It should be noted that due to difficulties that need to 
be faced in numerical research on condensing flows 
many correction factors modifying quantity J can be 
found in literature. However, their physical justification 
is very poor. One such coefficient is coefficient 8, which 
may limit the nucleation rate value. In practice, it is 
adopted based on experimental testing in convergent- 
divergent nozzles. The ranges of the selection of correc- 
tion factors for flows through nozzles are defined rela- 
tively well. For other flows through more sophisticated 
channels the situation is more complex due to differences 
in the nature of the expansion processes and a more lim- 
ited reference and experimental material. No correction 
coefficients are used in the numerical model applied 

herein, and coefficient £ in Equation (19) is equal to /. 


Droplet growth model 


Describing the processes of exchange between the 
droplet and the gaseous phase surrounding it, the basic 
principles of the molecular-kinetic mechanism of the 
interaction of the gaseous medium with the phase separa- 
tion surface need to be considered. The molecular-kinetic 
model is the case where the droplet free path length 1, is 
comparable to or bigger than the droplet diameter. Oth- 
erwise, the situation concerns a continuous medium 
model. The droplet free path-to-diameter ratio is called 
the Knudsen number expressed as follows: 


, K -n-R-T, 


s p 
Kn = = 21 
2er 2-r on) 
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In the case of water vapour condensation in moist air, 
the droplet free path value is much higher than the drop- 
let diameter: Kn>>1. The wetness mass fraction is here 
far smaller than in the case of pure water vapour. There- 
fore, the process of the heat and mass transfer between 
the gaseous and the liquid phase is based on the molecu- 
lar-kinetic theory [6]. Thus, the droplet growth equation 
may be written as: 


dr = a : Py — Ps (22) 


where a, is the condensation coefficient, here = 1, and 
Ps is the water vapour saturation pressure. 


Numerical Results 


Validation calculations - Condensation model valida- 
tion for moist air expansion in the Barltma nozzle 

The experimental testing of the moist air transonic 
condensing flow carried out by Dohrmann [3] in the 
Barltma hyperbolic nozzle is well known. The nozzle 
profile is described by the following equation: 


(23) 


where the nozzle height in the critical cross-section 
and the nozzle radius are y*=0.06 m and R*=0.2 m, re- 
spectively. 

For the purpose of the calculations presented herein, 
the flow with the inlet air relative moisture value ®ọ = 
29.4% at the pressure of the surroundings of 109 kPa and 
temperature of 300.6 K, was selected. The static pressure 
value at the level of 10 kPa is assumed at the outlet, 
which ensures a supersonic outflow. The numerical mesh 
adopted for the calculations is a structural, single-block 
mesh with the dimensions of 161<101*5. Fig. 1 presents 
the nozzle geometry and the distributions of the Mach 
number and the wetness mass fraction; Fig. 2 shows the 
distributions of static pressure and the wetness mass frac- 
tion along the nozzle central part. 

As it can be seen in Fig. 2, the comparison between 
computational and experimental results is satisfactory. 
Both the computational spontaneous condensation loca- 
tion and the condensation wave intensity are very close 
to experimental data. 

The moist air flow validation calculations were also 
performed for many other cases, both in nozzles and 
around blade profiles. They proved the developed in- 
house CFD code good capacity for modelling transonic 
condensing flows. 


Validation of the transonic flow with a normal shock 
wave 


The next validation test presented herein was con- 
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ducted using an experiment carried out by Sajben and 
Kroutil [11]. The test is very common and concerns the 
transonic flow through a de Laval “half-nozzle”. It was 
used to investigate the impact of back pressure in the 
nozzle on the normal shock wave location. Correct mod- 
elling of the shock wave location depends largely on the 
correct discretization of the viscous terms in the RANS 
equations and on the applied turbulence model. The com- 
putations were performed on a structural numerical mesh 
(Fig. 3). The boundary conditions adopted for the calcu- 
lations correspond to the shock wave location in the di- 
vergent part close to the critical cross-section y*. The 
pressure and the total temperature at the nozzle inlet were 
134470 Pa and 277.8 K, respectively, whereas static 
pressure at the outlet was 110660 Pa. 

The Mach number distribution in the nozzle is shown 
in Fig. 4, where the normal shock wave location, which 


Mach number contour lines (Ma = 0 + 2 with AMa = 0.2) 


Wetness contour lines (y = 0 + 7 gl/kg with Ay = 0.5 gl/kg) 
Fig. 1 Isolines of the Mach number and the mass moisture 
degree in the Barltma nozzle 


Calculated pressure 
Calculated wetness, y 
Measured pressure 


Fig. 2 Distributions of calculated and measured static pres- 
sure and of calculated wetness mass fraction along the 
nozzle central part 
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Mm 


-0.1 0 0.1 0.2 0.3 
x, m 


Fig. 3 Sajben's nozzle geometry and numerical mesh (3 


blocks: 31x151x3, 201x151x3 and 31x151x3). 


Fig. 4 Mach number contour lines (Ma=0+1.2 with AMa= 
0.05) for Sajben's nozzle 


for the adopted pressure ratio is close to the critical cross- 
section, can be seen clearly. 

The experiment provides a relatively large body of 
data that can be compared to the calculation results, such 
as the static pressure distribution along the nozzle upper 
and bottom walls, the velocity distribution for a few 
cross-sections of the nozzle, and the separation zone lo- 


Calculated pressure 
Measured pressure 


Bottom wall 


3 4 5 6 7 


F 2 at 8 i 2 
x/y* 
Upper wall 


Fig. 5 Distribution of calculated and measured static pressure 
on the nozzle bottom and upper wall 


Fig. 6 
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cation on the upper wall. With so much experimental data, 
the numerical code and the selected turbulence model 
may be validated reliably. Further on, Fig. 5 and Fig. 6 
show a comparison between the calculation results and 
the experimental data available from literature [11]. It 
can clearly be seen that the applied numerical code is 
able to model the transonic flow with a normal shock 
wave in Laval nozzles correctly. 


Calculations of the liquid phase evaporation on the 


normal shock wave at Barltm4 convergent-divergent 
nozzle 


For the Barltmaé CD nozzle, the outlet static pressure 
was raised to 70 kPa to generate a normal shock wave at 
the nozzle outlet. The generated shock wave was strong 
due to the fact that the expansion rate in the CD nozzle 
under consideration is high. There is an almost triple in- 
crease in pressure on the shock wave, which results in 
rapid and complete evaporation of the liquid phase cre- 
ated due to homogeneous condensation (Fig. 7) even if 
the relative moisture values at the nozzle inlet are high. 


Calculated axial velocity 
Measured axial velocity 
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It is known that on the normal shock wave there is a 
change in all parameters except total temperature (and for 
a perfect gas — also total enthalpy). In the case of moist 
air expansion, total temperature on the condensation 
wave increases initially, 47) con, due to the process of 
nucleation and droplets growth, whereas on the shock 
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80000 —--- Wetness, y fi 
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wave — there is a drop in total temperature, A To eva due to 
evaporation, and the total temperature value returns to the 
value at the nozzle inlet (Fig. 8). 

Fig. 9 and Fig. 10 present distributions of the Mach 
number and of the wetness mass fraction in a de Laval 
nozzle with a normal shock wave at the outlet for two 


18 0.015 
80000 
15 40.0125 
12 0.01 
60000 
i g 1 
o 0075 `- 
Sopi 0.0075 J 
0.6 + 40000 0.005 
0.3 0.0025 
20000 
o o1 0 0.1 02 
xm 
Dy = 70% 


Fig. 7 Distributions of static pressure, Mach number and wetness mass fraction along the nozzle central part 


Total temperature 
---- Static temperature 


xm 


Dy = 70% 


Fig. 8 Distributions of static and total temperature along the nozzle central part 


B, = 70% 


Fig. 9 Mach number contour lines (Ma = 0 + 1.5 with Ama = 0.1) 


Ð, = 70% 


Fig. 10 Moisture contour lines (y = 0 + 10g)/kg with Ay = 1 g)/kg) 
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values of relative moisture at the inlet. A clear impact can 
be seen of the relative moisture value at the inlet on the 
shock wave location — as the relative moisture value at 
the inlet rises, at other parameters at the nozzle inlet and 
outlet maintained at the same level, the shock wave is 
shifted downstream. Also the shock wave intensity de- 
creases as the relative moisture value gets higher. 


Sajben's transonic diffuser 


Similar analyses were conducted for the flow through 
Sajben's nozzle, for two values of relative moisture at the 
inlet and for the same as before values of total parameters 
at the inlet and of static pressure at the outlet. In this case, 
because the nozzle opening degree is smaller, the expan- 
sion rate is much slower and, consequently, the shock 
wave is much weaker (Fig. 11). Therefore, complete 
evaporation of the liquid phase occurs if the relative 
moisture value at the inlet is small, whereas for the rela- 
tive moisture value of 70% only a small part of the liquid 
phase evaporates on the shock wave. This is also the ef- 
fect of the shock wave weakening due to a rise in the 
moist air relative moisture. 

Fig. 12 and Fig. 13, respectively, show distributions of 
the Mach number and the mass moisture degree for the 
analysed cases of the moist air flow in Sajben's nozzle. 
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Fig. 11 Distributions of static pressure, the Mach number and 
the wetness mass fraction along Sajben's nozzle cen- 
tral part 
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Fig. 12 Mach number contour lines (Ma = 0 + 1.15 with Ama 
= 0.05) for Sajben's nozzle 


Fig. 13 Moisture contour lines (y = 0 + 22g/kg with Ay = 1 
g/kg) for Sajben's nozzle 


Conclusion 


This work considers the problem of the liquid phase 
evaporation due to condensation of water vapour con- 
tained in atmospheric air. The condensation process oc- 
curs as a result of moist air expansion in conver- 
gent-divergent nozzles. The capacity of the developed 
in-house CFD code for modelling this type of physical 
phenomena was tested. The following conclusions can be 
drawn based on the calculations: 

e The applied in-house CFD code models the water 
vapour condensing flow in convergent-divergent 
nozzles with a normal shock wave correctly. 

e In a Laval nozzle, complete evaporation of the liq- 
uid phase created due to homogeneous condensa- 
tion occurs on the shock wave almost in every case. 
Another water vapour condensation process may 
take place for example as a result of further ex- 
pansion of moist air. 

e Gentle and gradual evaporation of the liquid phase 
is possible in the case of weak shock waves. 

e A rise in the air relative humidity has an impact on 
the normal shock wave location and intensity — the 
wave is weakened and shifted downstream. 

e A rise in total temperature occurs on the condensa- 
tion wave, whereas on the shock wave, due to the 
liquid phase evaporation, the total temperature re- 
turns to the value that was initially at the nozzle 
inlet. 

Further works will focus on the developed methodol- 

ogy application for more complex geometries, as well as 
for flows with oblique shock waves. 
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